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Executive  Summary 


This  final  report  documents  the  major  developments  and  findings  during  the  grant  period  from 
March  2009  to  November  2012.  The  main  objective  of  this  project  was  to  develop  a  new 
discontinuous  formulation  named  correction  procedure  via  reconstruction  (CPR)  for  hyperbolic 
conservation  laws,  and  demonstrate  its  capability  for  the  Euler  and  Navier-stokes  equations  on 
hybrid  3D  unstructured  prismatic  and  tetrahedral  grids.  The  CPR  method  can  unify  several 
popular  high  order  methods  including  the  discontinuous  Galerkin,  the  spectral  volume  and 
spectral  difference  methods  into  an  efficient  differential  form  without  explicit  volume  or  surface 
integrals.  By  selecting  the  solution  points  to  coincide  with  the  flux  points,  solution  reconstruction 
can  be  completely  avoided. 

We  successfully  fulfilled  the  main  objective  in  the  present  study.  More  specifically,  we  achieved 
the  following  accomplishments: 

•  Extended  the  CPR  formulation  to  3D  hybrid  meshes,  including  tetrahedral,  hexahedral, 
prismatic  elements; 

•  Extended  the  CPR  formulation  to  the  Navier-Stokes  equations  on  hybrid  elements,  and 
demonstrate  the  method  for  benchmark  3D  problems 

•  Implemented  the  CPR  method  on  clusters  of  CPUs  and  GPUs,  and  achieved  up  to  two  orders 
of  magnitude  speedup  on  the  GPU  than  the  CPU 

•  Extended  the  method  for  dynamic  moving  grids  satisfying  the  so-called  geometric 
conservation  laws,  and  demonstrated  the  capability  for  bio-inspired  flow  problems 

•  Implemented  solution-based  hp-adaptations  using  a  variety  of  adaptation  criteria  including 
residual,  adjoint  and  entropy  based  adaptation  criteria. 

In  the  First  International  Workshop  on  High-Order  CFD  Methods,  it  was  demonstrated  through 
many  benchmark  cases  that  high-order  methods  are  capable  of  achieving  the  same  error  with  less 
CPU  time  than  low  order  methods. 

The  grant  supported  1  postdoc  research  associate,  Dr.  T.  Haga,  and  several  PhD  students,  H. 
Gao,  Y.  Zhou  and  L.  Shi.  The  publications  partly  supported  by  the  grant  are  listed  below: 

•  Z.J.  Wang,  K.J.  Fidkowski,  R.  Abgrall,  F.  Bassi,  D.  Caraeni,  A.  Cary,  H.  Deconinck,  R. 
Hartmann,  K.  Hillewaert,  H.T.  Huynh,  N.  Kroll,  G.  May,  P-O.  Persson,  B.  van  Leer,  and  M. 
Visbal.  “High-Order  CFD  Methods:  Current  Status  and  Perspective,”  International  Journal 
for  Numerical  Methods  in  Fluids,  2012,  Accepted. 

•  H.  Gao  and  Z.J.  Wang,  “A  Conservative  Correction  Procedure  via  Reconstruction 
Formulation  with  the  Chain-Rule  Divergence  Evaluation”,  J.  Computational  Physics  232,  7- 
13  (2013). 

•  Y.  Li  and  Z.J.  Wang,  “An  Optimized  Correction  Procedure  via  Reconstruction  Formulation 
for  Broadband  Wave  Computation”,  Communications  in  Computational  Physics,  Vol.  13, 
No.  5,  pp.  1265-1291  (2013). 

•  M.  Yu,  Z.  J.  Wang  and  H.  Hu,  “Airfoil  Thickness  Effects  on  the  Thrust  Generation  of 
Plunging  Airfoils”,  J.  of  Aircraft,  accepted. 


2 


•  H.  Gao,  Z.J.  Wang  and  H.T.  Huynh,  “Differential  Formulation  of  Discontinuous  Galerkin 
and  Related  Methods  for  the  Navier-Stokes  Equations”,  Communications  in  Computational 
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1.  Introduction 


Advantages  of  high-order  methods  are  well  recognized  in  the  computational  fluid  dynamics 
(CFD)  community  especially  for  aeroacoustic  noise  predictions,  vortex  dominated  flows,  large 
eddy  simulation  and  direct  numerical  simulation  (DNS)  of  turbulent  flows.  Since  the  truncation 
error  of  a  high-order  method  decreases  more  rapidly  than  that  of  a  lower  order  method,  the  more 
stringent  the  accuracy  requirement  is,  the  more  efficient  a  high-order  method  becomes  in 
computational  cost.  For  the  practical  use  in  industries,  lower  order  (1st  or  2nd)  unstructured  grid 
methods  are  usually  employed  for  the  reason  of  superior  geometrical  flexibility  and  robustness. 
However,  these  methods  are  likely  too  dissipative  to  capture  small  vortex  structures  in  turbulent 
flows  and  are  often  not  capable  of  obtaining  grid  converged  solutions.  Increased  prediction 
accuracy  is  often  required  for  many  aerodynamic  problems  with  both  complex  physics  and 
geometry,  such  as  helicopter  blade  vortex  interactions,  flow  over  high  lift  devices,  and  aero¬ 
acoustic  noise  generated  by  the  landing  gear. 

In  the  past  decades,  there  has  been  significant  progress  in  developing  high-order  methods 
capable  of  solving  the  Navier-Stokes  (NS)  equations  on  unstructured  grids.  For  compressible 
flow  computations  in  aerospace  applications,  the  discontinuous  Galerkin  (DG)  method 
[31,4,5,7,1,2,29,42]  has  attracted  intensive  interest.  One  particular  feature  of  the  DG  method  is 
the  discontinuous  solution  space  of  high-order  approximations  for  each  element,  which  allows 
the  scheme  to  be  very  flexible  in  dealing  with  complex  configuration  and  in  accommodating 
solution  based  adaptations.  Other  methods  assuming  element-wise  discontinuous  solution  are 
staggered-grid  (SG)  multi-domain  spectral  method  [21],  spectral  volume  (SV)  [43,46- 
48,25,39,14,12]  and  spectral  difference  (SD)  [23,24,28,40]  methods.  Another  notable  feature  that 
is  common  among  these  methods  is  the  use  of  one  of  the  Riemann  solvers  [33,32,30,19,22]  to 
compute  unique  fluxes  at  element  interfaces  to  incorporate  “upwinding”  characteristics  of  wave 
propagation,  similar  to  the  Godunov  type  finite  volume  method  [11,41],  The  main  difference 
among  these  methods  lies  in  how  the  governing  equations  are  discretized  and  the  degrees-of- 
freedom  (DOFs)  are  chosen.  The  DG  method  is  based  on  the  weighted  residual  form.  Various 
types  of  DG  schemes  are  derived  with  different  choice  of  DOFs.  Depending  on  how  the  DOFs 
are  defined,  DG  schemes  can  be  further  divided  into  modal  and  nodal  approaches.  The  SV 
method  is  discretized  in  the  integral  form  similar  to  the  finite  volume  method  and  the  DOFs  are 
always  the  sub-cell  or  control  volume  (CV)  averages.  The  SG/SD  method  is  based  on  the 
differential  form  without  any  integration  and  the  DOFs  are  chosen  as  the  nodal  values  within 
each  element.  More  comprehensive  reviews  of  these  methods  are  given  in  [44] 

Recently,  a  novel  formulation  named  CPR  (correction  procedure  via  reconstruction)  was 
developed  by  Huynh  [17,18]  for  ID  conservation  laws,  and  extended  to  simplex  and  hybrid 
meshes  by  Wang  and  Gao  [45].  The  CPR  method  is  based  on  a  nodal  differential  form,  with  an 
element-wise  discontinuous  polynomial  solution  space.  The  solution  polynomial  is  interpolated 
from  the  solutions  at  a  set  of  solution  points.  This  formulation  has  some  remarkable  properties. 
The  framework  is  easy  to  understand,  efficient  to  implement  and  recovers  several  known 
methods  such  as  the  DG,  SG  or  the  SV/SD  methods.  Furthermore,  by  choosing  the  solution 
points  to  coincide  with  the  flux  points,  the  reconstruction  of  solution  polynomials  to  calculate  the 
residual  can  be  completely  avoided.  The  DG  scheme  derived  through  the  CPR  framework  is 
probably  the  simplest  and  most  efficient  amongst  all  DG  formulations  since  explicit  integrations 
are  avoided.  In  a  recent  study  [9],  the  CPR  method  has  been  extended  to  the  Navier-Stokes 
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equations  on  2D  mixed  meshes.  These  successful  developments  laid  a  solid  foundation  for  its 
efficient  implementation  and  demonstration  on  arbitrary  grids  in  3D. 

Hybrid  elements  such  as  hex,  prism,  pyramid  and  tetrahedron  will  provide  great  geometrical 
flexibility  for  practical  problems  in  3D.  In  particular,  for  high  Reynolds  number  flows  in 
aerodynamic  applications,  prismatic  cells  have  the  advantages  in  accuracy  and  computational 
cost  to  resolve  boundary  layers  near  the  wall.  There  have  been  several  attempts  to  develop  the 
DG  method  on  arbitrary  grid  elements.  In  [34]  different  types  of  elements  such  as  hex,  prism  and 
pyramid  are  projected  onto  a  reference  cube  using  collapsed  Cartesian  coordinates  and 
hierarchical  basis  functions  over  the  cube  are  used.  Luo  et.al.  [26]  presented  a  different  approach 
based  on  the  Taylor  series  expansion  at  the  center  of  arbitrary  element.  Gassner  et.al.  [10]  used 
polymorphic  nodal  element  in  the  modal  based  formulation  to  reduce  the  cost  of  numerical 
integrations.  However  one  obvious  shortcoming  of  these  formulations  is  the  high  computational 
cost  of  the  surface  and  volume  integrations  coming  from  the  weighted  residual  formulation. 
Another  difficulty  is  the  treatment  of  curvilinear  boundary  elements.  The  simple  formulation  of 
the  CPR  method  is  expected  to  alleviate  the  computational  costs  and  facilitate  the  treatment  of 
curved  wall  surfaces  in  a  straightforward  fashion. 

In  the  present  study,  we  develop  the  CPR  for  solving  the  Euler  and  Navier-Stokes  equations 
on  3D  mixed  meshes.  For  the  current  implementation,  tetrahedral  and  prismatic  elements  are 
considered  with  the  intention  to  resolve  viscous  boundary  layer  flows  efficiently.  The  remainder 
of  this  article  is  organized  as  follows.  The  basic  formulation  of  the  CPR  method  is  described  in 
the  next  section.  In  section  3,  The  discretization  of  the  compressible  Navier-Stokes  equations  is 
derived  in  the  CPR  framework.  Subsequently,  we  discuss  how  to  implement  the  CPR  method 
efficiently  in  each  particular  element  with  curvilinear  geometry  in  section  4.  Section  5  presents 
the  computational  results  for  several  benchmark  problems,  including  accuracy  studies  on  mixed 
unstructured  grids.  Conclusions  for  the  present  study  and  possible  future  works  are  summarize  in 
section  6. 

2.  Review  of  the  Correction  Procedure  via  Reconstruction  Formulation 


We  first  review  the  CPR  formulation  for  a  hyperbolic  conservation  law,  which  can  be  written 
as 


^  +  V*F(0  =  O,  (2.1) 

dt 

with  suitable  initial  and  boundary  conditions.  Q  is  the  vector  of  conserved  variables,  and  F  is  the 
flux  vector.  Assume  that  the  computational  domain  is  discretized  into  N  non-overlapping 
elements  {VJ.  The  weighted  residual  form  of  (2.1)  on  element  vi  can  be  derived  by  multiplying 
(2.1)  by  an  arbitrary  weighting  or  test  function  W  and  integrating  over  v;, 

j  ^-WdV  +  j  WF(Q )  •  MS  -  j  Vlk  •  F(Q)dW  =  0.  (2.2) 

v i  av,  v, 

Let  Q1'  be  an  approximate  solution  to  Q  at  element  Vr  We  assume  that  the  solution  belongs  to  the 
space  of  polynomials  of  degree  k  or  less,  i.e.,  Qht  e  Pk,  within  each  element  without  continuity 
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requirement  across  element  interfaces.  Then,  we  require  that  the  numerical  solution  Q1.'  must 
satisfy  (2.2),  i.e.. 


j  - 'Q-WdV+  J  WF(Q.  )»ndS-  j  VWF(Q?)dW  =0.  (2.3) 

v,  V, 

Because  the  approximated  solution  is  in  general  discontinuous  across  element  interfaces,  the 
fluxes  at  the  interfaces  are  not  well  defined.  To  evaluate  a  unique  flux  and  also  to  provide 
element  coupling,  a  common  Riemann  flux  is  used  to  replace  the  normal  flux,  i.e., 

Fn{Qhi)  =  F{Ql;)*n~F\Qhi,Ql,n),  (2.4) 

where  Q'‘+  is  the  solution  from  outside  of  the  current  element  Vr  Thus,  Eq.  (2.3)  becomes 

J  ^Q-WdV+  j  WFn(Q';,Q';+,n)dS-  J  VI VF(Qk)dW  =0.  (2.5) 

V,  dVi  V, 

If  the  space  of  W  is  chosen  to  be  the  same  as  the  solution  space,  Eq.  (2.5)  is  equivalent  to  the  DG 
formulation.  For  the  sake  of  a  much  simpler  formulation,  we  wish  to  eliminate  the  test  function. 
Applying  integration  by  parts  to  the  last  term  of  (2.5),  we  obtain 

J  ^Q-WdV  +  J  WV  •  F(Q.‘)dW  +  J  W  [F"(Q'' -  F"(gf  )]d£  =  0.  (2.6) 

v i  dt  v,  av, 

Note  that  the  surface  integral  coming  from  the  last  term  of  (2.5)  is  evaluated  by  only  using  the 
internal  solution.  The  last  term  of  (2.6)  can  be  viewed  as  a  penalty  term,  i.e.,  penalizing  the 
normal  flux  differences  [F]  =  Fn  (Qk  ,Qk+,n)  -  Fn  (Qk) .  Let  us  introduce  a  “correction  field”  S;  e  Pk , 
which  is  determined  from  the  following  relation  defining  the  so-called  “lifting  operator”  for  [F]. 

J  WStdV  =  J  W[F\dS.  (2.7) 

v,  av, 


Substituting  (2.7)  into  (2.6),  we  obtain 


J 


+  V»F(Q’il)  +  SiWdV=  0. 


(3.8) 


In  the  present  study,  in  order  to  simplify  the  derivation  we  also  approximate  the  flux  divergence 
by  polynomials  of  degree  k  or  less,  i.e.  V  •  F( Q. )  e  Pk.  If  W  is  selected  such  that  a  unique  solution 
exists,  (3.8)  is  equivalent  to 


dt 


+  V*F(<2,")  +  3  =  0, 


(2.9) 
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i.e.,  (2.9)  is  satisfied  everywhere  in  element  Vr  With  the  definition  of  a  correction  field,  we  have 
successfully  reduced  the  weighted  residual  formulation  to  an  equivalent  simple  differential  form, 
which  does  not  involve  any  explicit  surface  or  volume  integrals. 

To  find  the  approximate  solution  gf,  let  the  DOFs  be  the  solution  values  at  a  set  of  points 
{rtj },  named  solution  points  (SPs).  Then  equation  (2.9)  must  hold  at  the  SPs,  i.e., 


Ml 

dt 


+  V*F(Q';I)  +  SI]=Q. 


(2.10) 


Obviously  the  correction  field  5i  can  be  expressed  in  terms  of  8i .  using  a  Lagrange  interpolation 
on  the  SPs,  i.e., 


4  =  I>5 


(2.11) 


where  LSP  e  Pk  are  the  Lagrange  polynomials  based  on  the  SPs.  In  the  case  of  a  non-linear  flux 
vector,  F(Q.)  are  not  polynomials  in  general.  In  the  present  study,  we  approximate  F  (£)'')  by 
polynomials  of  degree  k  to  evaluate  RHS  of  eq.  (2.7).  Thereby,  we  assume  that  the  flux 
difference  [F]  is  a  polynomial  on  the  faces  of  the  element,  and  can  be  determined  based  on 
values  of  [F]  at  a  set  of  flux  points  (FPs)  { rf , )  using  a  Lagrange  interpolation,  i.e., 


[^]/=Z  L7(rfJ)[F]fJ,  (2.12) 

/ 

where  Lw  e  Pk  are  the  Lagrange  polynomials  based  on  the  FPs.  Then,  if  the  locations  of  the 
solution  and  flux  points  are  specified  and  the  weighting  function  W  is  chosen  so  as  to  have  the 
same  dimension  as  the  correction  field  Sp  Stj  can  be  uniquely  defined  by  solving  the  linear 

system  derived  from  eq.  (2.7).  For  simplex  elements  with  straight  faces,  it  can  be  expressed  in 
the  following  formula 


4,  - 


Z  Z<W^/A 

/s3Vj  l 


(2.13) 


where  a  are  constant  coefficients  independent  of  the  solution.  Substituting  (2.13)  into  (2.10) 
we  obtain  the  following  equation 


M 


ajfl[F]flSf=  0. 


(2.14) 


One  can  clearly  see  that  this  is  a  collocation-like  formulation  with  penalty-like  term  that  comes 
from  the  element-wise  correction  polynomial  to  provide  the  coupling  between  elements.  It  can  be 
shown  that  the  location  of  SPs  does  not  affect  the  numerical  scheme  for  linear  conservation  laws 
[40,17].  For  efficiency,  the  solution  points  are  always  chosen  to  coincide  with  the  flux  points. 
Therefore,  any  data  interpolation  is  no  longer  needed  for  flux  calculation,  which  dramatically 
reduces  the  computational  cost.  Any  convergent  nodal  sets  with  enough  points  at  the  element 
interface  are  good  candidates,  e.g.,  those  found  in  [3,15,49]. 
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Finally  we  want  to  make  a  remark  on  the  relationship  between  the  CPR  formulation  and  other 
methods  including  DG,  SV  and  SD  methods.  Starting  from  the  weighted  residual  form  of  the 
governing  equations,  different  formulations  can  be  derived  depending  on  the  weighting  function. 
For  example,  a  nodal  DG  formulation  is  obtained  by  choosing  weighting  functions  to  be 
Lagrange  polynomials,  and  a  SV  formulation  is  obtained  by  defining  weighting  functions  as 
piecewise  constant  at  the  sub-cells.  As  a  result,  the  only  difference  between  those  schemes 
appears  in  the  correction  coefficients.  In  the  original  work  [45],  it  was  shown  that  the  resulting 
CPR  scheme  is  basically  conservative  by  using  the  correction  coefficients  for  the  DG,  SV  and 
SD  scheme.  In  this  study,  we  choose  the  weighting  function  to  be  one  of  the  Lagrange 
polynomials  based  on  the  SPs,  i.e.,  eq.  (2.9)  is  identical  to  the  DG  formulation. 

3.  Discretization  of  the  Navier-Stokes  Equations 


3.1.  Governing  Equations 

The  3D  compressible  Navier-Stokes  equations  can  be  written  as  a  system  of  partial 
differential  equations  in  conservation  form: 

^  +  V  •  (Fc  (£>)  -  Fv  (Q,  V  £>))=  0,  (3.1) 


where  Q  ,  Fc  =  [Fcx,Fcy,Fcz]  and  Fv  =  [FX,FJ,FZ]  denote  the  conservative  state  vector,  the  inviscid  and 
the  viscous  flux  vectors,  respectively,  and  are  given  by 
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K(e  +  p)wj 

0  ] 

°  ] 

'  0 

L* 

Lx 

.  Ky  = 

Tyy 

.  K  = 

Tzy 

Txz 

Tyz 

Tzz 

Ku^xx  +  vrx7  +  wZxz-qx) 

yUTyx  +  VTyy  +  WTyz-qy  ) 

(3.2) 


(3.3) 


where  p  is  the  density,  v  =  [u,v,w]  are  the  velocity  vector,  p  is  the  pressure,  e  is  the  total  energy 
per  unit  volume.  The  viscous  stress  tensor  can  be  represented  as 


T  = 


+  (Vv)  --(V»LK 


(3.4) 


where  p  is  the  molecular  viscosity  coefficient,  /is  the  unit  tensor.  The  heat  flux  is  given  as 


q  =—c—VT 

p  Pr 


(3.5) 
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Here,  Cp  is  the  specific  heat  capacity  at  constant  pressure  and  T  is  the  temperature.  The  Prandtl 

number  Pr  is  assumed  to  be  a  constant  of  0.72  in  this  study.  For  a  perfect  gas,  the  pressure  is 
related  to  the  total  energy  e  by 


e  =  — —  +  —  p(u2  +  v2  +  w2\  (3.6) 

y- 1  2  v  7 

The  specific  heat  ratio  y  is  set  to  be  a  constant,  1.4  for  air.  The  computations  for  solving  the 
Euler  equations  are  performed  by  omitting  the  viscous  flux. 

3.2.  CPR  Formulation  of  the  Navier-Stokes  Equations 

In  order  to  discretize  the  Navier-Stokes  equations,  we  follow  a  mixed  formulation  that  is 
commonly  used  for  the  DG  method[2,  6].  By  introducing  a  new  variable  R  =  VQ,  Eq.  (3.1)  is 
rewritten  in  a  first  order  system  as 


^  +  V(Fc(Q)-Fv(Q,R))=0, 


(3.7) 


R  =  VQ. 


(3.8) 


According  to  the  CPR  formulation  by  assuming  Q.,R*  e  Pk  on  discretized  elements  \V:},  we 
obtain 


^fi  +  v*(f.<a‘,)-f.(e‘,.*‘,))+4j  I  <3'9) 


HfedVi  l 


<-<vcVp  2 


i I  f&dVi  l 


where  [FJ  =  Fc”(Q?,Q?+,n)  -  Fcn(Q»), 

[Q]  =  Q(Q'\  Q'l  )-£>■■ 


[Fv\ -  F^QlQlyQlVQln) -  F"(Q‘,R’) 


(3.10) 

and 


3.2.1.  Inviscid  Flux  Calculation 

Here  we  consider  the  inviscid  flux.  We  need  to  discretize  the  internal  flux  divergence  and  the 
common  flux  at  the  interface  in  (3.9).  Instead  of  approximating  the  inviscid  flux  by  the  Lagrange 
interpolation  on  the  SPs,  the  flux  divergence  is  calculated  “exactly”  at  the  solution  points  with 
the  chain  rule  (CR)  approach 


V‘C(<2M  = 


dFMj) 

a<2 


'VO" 


(3.11) 
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■Jjn 

where  is  the  inviscid  flux  Jacobian  matrix.  Note  that  V  •  FJQ'.  )  is  generally  not  a  degree  k 
dQ 

polynomial,  but  it  can  be  approximated  by  the  Lagrange  polynomial  of  degree  k  using  the  flux 
divergence  at  the  solution  points,  i.e., 

v  •  Fc(Q>‘)  «  V  •  FcCR(Gf)  =  Ojf.)V  •  Fe(Q*j),  (3.12) 

j 

This  implies  that  the  flux  vector  FcCR  ( Q. )  belongs  to  Pk+I  which  is  one  degree  higher  than  the 
approximation  of  F(Qht)  e  Pk  used  in  the  correction  term.  The  chain  rule  approach  is  known  to  be 
more  accurate,  though  the  resulting  scheme  is  not  strictly  conservative  due  to  the  inconsistency 
between  the  approximated  flux  vectors  in  the  flux  divergence  term  and  the  correction  term. 
However,  it  is  also  known  that  the  mass  conservation  error  is  extremely  small  [45].  If 
conservation  is  absolutely  required,  one  can  use  the  Lagrange  polynomial  approach  [17]. 

The  common  inviscid  flux  in  [FJ  can  be  obtained  with  any  Riemann  solver.  In  this  paper,  the 
Roe  flux  [32]  is  used  for  all  the  cases. 


3.2.2.  Viscous  Flux  Calculation 

In  the  present  study,  we  employ  the  BR2  scheme  [2]  to  discretize  the  viscous  flux.  In  (3.10), 
the  common  solution  Qf ;  is  simply  the  average  of  the  solutions  at  both  sides  of  /.  The  viscous 

fluxes  at  the  solution  points  are  evaluated  by  Fv(<2*. ,/?''.).  Then  the  viscous  flux  divergence  is 

obtained  through  the  Lagrange  interpolation  rather  than  the  CR  approach.  In  the  correction  term, 
the  common  viscous  flux  also  needs  to  be  determined.  Besides  the  common  solution,  we  also 
need  to  define  a  common  gradient  on  face/.  The  common  gradient  is  evaluated  as 

v2"r  =  j  (VC;.,  +  r~j  +  V0* ,  +  rjj)  (3.13) 


where  VQfl  and  V<2/  are  the  gradients  of  the  solution  from  the  left  and  right  cells,  while  rfl  and 
r/  are  the  local  lifting  correction  to  the  gradients  only  due  to  the  common  solution  on  face/ 


Iv1 


teiy, ,„(+«/>/> 


(3.14) 


where  m  is  the  index  for  the  flux  points  on /and  hf  is  the  unit  normal  vector  directing  from  left 

to  right.  Note  that  there  is  no  summation  over  all  faces  of  the  element  in  eq.  (3.14)  in  order  to 
assure  local  property  of  the  BR2  scheme. 


4.  Discretization  on  Mixed  Grids  with  Curved  Boundary 


It  can  be  observed  that  (2.9)  is  valid  for  arbitrary  types  of  elements  besides  triangles  and 
tetrahedrons.  The  current  development  for  3D  hybrid  meshes  accommodates  two  kinds  of 
element  shapes,  i.e.,  tetrahedron  and  triangular  prism.  Other  types  of  element  such  as  hexahedron 
and  pyramid  will  be  developed  in  the  near  future.  The  use  of  prismatic  cells  in  addition  to 
tetrahedral  cells  has  the  advantages  in  both  accuracy  and  computational  costs  to  resolve 
boundary  layers  near  solid  walls.  In  order  to  achieve  an  efficient  implementation,  all  elements 
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are  transformed  from  the  physical  domain  (x,  y,  z)  into  a  corresponding  standard  element  in  the 
computational  domain  (<£  tj,  £)  as  shown  in  Fig.  1.  Here  we  consider  the  transformations  for  the 
elements  with  curved  sides  (faces  and  edges).  The  discretization  for  the  curved  elements  is 
conducted  in  the  same  way  as  the  straight  sided  elements  by  applying  the  CPR  formulation  in  the 
standard  elements.  In  the  present  study,  a  quadratic  triangular  face  is  employed  to  represent 
curved  wall  boundaries.  For  the  sake  of  computational  efficiency,  the  quadratic  representation  is 
adopted  for  only  one  of  the  faces  of  tetrahedron  which  will  be  attached  to  the  wall  in  inviscid 
flows,  and  for  only  two  triangular  faces  of  prism  which  will  be  used  in  the  thin  layers  of  prism 
cells  to  assure  the  quality  of  the  element  shape  especially  in  high  Reynolds  number  flows. 

Based  on  a  set  of  locations  of  nodes  defining  the  shape  of  element,  a  set  of  shape  functions 
can  be  obtained  [50].  Once  the  shape  functions  are  given,  the  transformation  can  be 

written  as 


(  \ 
X 

r  1 

(  \ 
xi 

y 

J 

II 

yt 

UJ 

*4  1 

UJ 

where  K  is  the  number  of  points  used  to  define  the  physical  element,  (jcf,y,.,z,.)  are  the  Cartesian 
coordinates  of  those  points.  For  the  transformation  given  in  (4.1),  the  Jacobian  matrix  J  takes  the 
following  form 


y # 


XV  XC 

y,  yc 

Zq  Zg  _ 


(4.2) 


For  a  non-singular  transformation,  its  inverse  transformation  must  also  exist,  and  the  Jacobian 
matrices  are  related  to  each  other  according  to 


4M 

d(x,y,z) 
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6 
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£ 

% 

U 


=j~l. 


(4.3) 


The  governing  equations  in  the  physical  domain  are  then  transformed  into  the  computational 
domain  (standard  element),  and  the  transformed  equations  take  the  following  form 


dQ  dFi  dF ”  dFc  n 
dt  dr]  dC 


(4.4) 


where 
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Figure  1.  Transformations  of  a  curve  boundary  tetrahedral  and  prismatic  cell  to  the  standard 

Q  =  \J\-Q, 

Ft=\j\.(txF>  +  ZyF’  +  tr)  (4.5) 

F”=\j\(rjxFx  +  JJyFy  +  JjzFz) 

F(  =  M  •  (CJ?X  +  CyFy  +  £ZFZ) 


Let  Stl  =  \j\(i]x,i)y,i]7)  and  Sc  =\j\(Cx,Cy,Q-  Then  we  have  F*  =F*Sr  F"  =F*S„  and 

Fc  =  F»Sr  In  our  implementation,  |/|,  Sr  Sr/  and  S(  are  stored  at  the  solution  points.  Note  that 

here  we  consider  the  Euler  equations  as  the  governing  equations  for  brevity’s  sake.  Extending 
the  following  discretization  to  the  Navier-Stokes  equations  is  straightforward. 

Discretization  on  a  Standard  Tetrahedron 

On  a  standard  tetrahedron,  the  CPR  formulation  in  eq.  (2.14)  can  be  rewritten  as 

+  1  X  =  0,  (4.6) 

M  \V  |  fedV  l 

where  superscript  (0  means  the  variables  or  operations  evaluated  on  the  computational  domain. 
For  example,  [F(f)]  are  the  normal  jumps  of  the  transformed  fluxes  across  the  faces  of  the 
standard  element.  The  transformed  normal  flux  can  be  expressed  in  terms  of  the  flux  in  the 
physical  space  as 


?(£)  - 
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where  n^]  =  [n^ ,nv ,n(]  is  a  unit  normal  vector  on  a  straight  face  of  the  standard  element,  and  Sn  is 
a  normal  vector  on  a  face  in  the  physical  space  defined  as  sn  =  +  St]nn  +  Scnc. 

Note  that  solving  Eq.  (4.4),  Q  =  \j\-Q  are  the  solution  unknowns,  and  are  assumed  to  be  degree 
k  polynomials  in  the  computational  domain  instead  of  Q.  As  a  result,  the  derivatives  of  Q  should 
be  calculated  in  the  following  way. 


dQ 


u\ 


4]\Q)  47l 


dQ 

drj 
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fV|e)_4 !\0 

drj  drj 


dQ 

dc 
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d(\AQ)  M 
dc  dC 


(4.8) 


In  3D,  to  construct  a  complete  polynomial  of  degree  k,  at  least  k(k+\)(k+2)/(  1x2x3)  SPs  need 
to  be  chosen.  In  order  to  achieve  the  most  efficient  implementation,  SPs  on  edges  are  chosen  to 
be  the  Legendre-Gauss  Lobatto  (LGL)  points.  For  4th-  or  higher  order  schemes,  nodes  inside  the 
boundary  triangle  are  chosen  from  [15].  For  5th-  or  higher  order  schemes,  nodes  inside  the 
tetrahedron  are  chosen  from  [49].  The  nodal  set  of  the  4th-order  CPR  scheme  is  shown  in  Fig.  2. 
Note  that  the  flux  difference  at  a  flux  point  corrects  all  solution  points  as  shown  in  (4.6). 


Discretization  on  Standard  Prism 

For  a  standard  triangular  prism,  the  solution  polynomial  can  be  expressed  as  a  tensor  product 
of  a  ID  and  2D  Lagrange  polynomial,  i.e., 

G*(£* 1,0 = YLokjM&WQ’  (4-9) 

k=  1  j=  1 


where  g*  are  the  state  variables  at  the  solution  point  (j,k),  with  j  the  index  in  ^-ij  plane  and  k 
the  index  in  C,  direction,  L^,rj)  is  a  2D  Lagrange  polynomial  in  a  triangle  and  Lk(£)  is  a  ID 

Lagrange  polynomial  in  a  segment.  Figure  3  shows  the  locations  of  the  solution  points  for  k= 3. 
The  nodal  sets  on  the  edge  and  the  triangle  are  chosen  in  the  same  manner  with  the  tetrahedral 
element. 

For  the  extension  of  the  CPR  method  to  2D  quadratic  elements  [17],  the  solution  procedure 
reduces  to  a  series  of  one-dimensional  operations,  i.e.,  the  solution  polynomial  is  represented  as  a 
tensor  product  of  ID  Lagrange  polynomials  and  the  correction  due  to  flux  difference  is 
performed  in  a  one-dimensional  manner.  The  CPR  formulation  for  a  standard  prism  can  be 
derived  in  an  analogous  fashion  as 


v( ^ 

VTri  fedVTri  l 


(O 

/ 


(4.10) 


+[Fc^j,7ij-l)-F^^j,rlj-l)yL  (Q+  [Ff(#y,^,l)-Ff(^,7y,l)>;  (Q  =  0. 


The  correction  process  is  done  in  a  decoupled  manner.  The  third  term  is  the  correction  of  the 
flux  components  in  £  and  ij  direction,  which  is  computed  on  a  plane  with  fixed  £  =  C,  This  is 
nothing  but  the  correction  used  in  the  2D  CPR  method  for  a  triangular  element.  In  eq.  (4.10),  Vt„ 
is  the  area  of  triangle,  5/  the  length  of  the  edge  /  and  /  the  index  for  flux  points  on  /.  Note  that, 
corrects  only  the  solution  points  on  the  triangle  with  fixed  k  instead  of  all  solution 

points  in  the  element  as  shown  in  Fig.  3  (a).  The  last  two  terms  denote  the  correction  in  the  g 
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Figure  2.  Solution  points  in  the  standard  tetrahedral  cell  for  degree  k=3  polynomial  (only 

points  on  the  visible  faces  are  shown). 

direction,  which  is  evaluated  by  the  ID  CPR  method  [17].  gi  and  gR  are  the  correction  functions 
for  the  left  and  right  end  points  of  a  segment.  The  flux  difference  at  an  end  point  corrects  only 
the  solution  points  on  the  segment  with  fixed  j  as  shown  in  Fig.  3  (b).  For  prism  cells,  the 
number  of  solution  points  corrected  by  a  flux  point  is  smaller  than  the  one  for  tetrahedral  cells 
due  to  the  decoupled  correction  procedure.  Hence,  the  method  for  prisms  is  more  efficient  per 
DOF  than  for  tetrahedrons.  This  decoupled  procedure  also  facilitates  the  implementation 
employing  different  degrees  of  polynomials  in  %-ij  and  ^directions  to  adapt  to  flow  features.  An 
attempt  employing  higher  order  polynomials  in  the  wall  normal  direction  to  resolve  the  boundary 
layer  with  coarser  prism  cells  is  shown  in  a  later  section. 

In  order  to  simplify  the  implementation  for  mixed  grids,  we  assume  the  polynomial  degree  k 
to  be  the  same  for  both  the  tetrahedral  and  prismatic  elements.  Furthermore,  the  flux  points  along 
the  element  interfaces  are  required  to  match  each  other.  In  the  present  implementation,  the  flux 
points  are  selected  to  be  the  LGL  points  at  each  edge  for  all  tetrahedral  and  prismatic  elements. 
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5.  Numerical  Results 

5.1.  Test  cases  for  the  Euler  Equations 

5.1.1.  Accuracy  Study  with  Vortex  Evolution  Problem 

To  assess  the  order  of  accuracy  of  the  developed  method,  the  propagation  of  an  isentropic 
vortex  in  inviscid  flow  is  computed  with  successive  grid  refinement.  This  is  an  idealized  problem 
for  the  Euler  equations  in  2D  used  by  Shu  [35],  Here  we  consider  simple  extension  of  this 
problem  to  the  3D  domain  [0,  10]x[0,  10]x[0,  10].  The  mean  flow  is  {p,  u,  v,  w,  p}={  1,  1,  1,  0, 
1}.  An  isotropic  vortex  is  then  added  to  the  mean  flow,  i.e.,  with  perturbations  in  u,  v,  and 
temperature  T  =  pip,  and  no  perturbation  in  entropy  S  =  p!  pr\ 


c  C 


Figure  3.  Solution  points  in  the  standard  prism  cell  for  degree  k= 3  polynomial  (only  points  on 
the  visible  faces  are  shown),  (a)  shows  the  correction  in  the  <f  and  //  derections.  (b)  shows  the 


Su  =  -yJ-e™**\  &  =  x—e™^\ 

2n  2  n 

(r-y  — F. 

8pr2 


Sw  =  0 


(5.1) 


where  72  =  I2  +  y2,  x  =  x-5,  y  =  y- 5 ,  and  the  vortex  strength  e=  5 .  If  the  computational  domain 
is  infinitely  big,  the  exact  solution  of  the  Euler  equations  with  the  above  initial  condition  is  just 
the  passive  convection  of  the  isentropic  vortex  with  the  mean  velocity  (1,1,  0).  In  the  numerical 
simulation,  we  impose  the  exact  solution  on  the  boundaries. 

The  computations  are  carried  out  until  t= 2  on  two  different  types  of  grids,  tetrahedral  meshes 
and  prismatic  meshes.  In  generating  computational  grids,  first  an  equidistant  Cartesian  grid  of 
NxNxN  cells  is  assumed  for  the  cubic  domain  and  each  cell  is  further  divided  into  six 
tetrahedrons  or  two  prisms.  Three  different  grids  are  employed  with  N=10,  20  and  40  for  each 
type  of  cell.  For  the  time  integration,  a  3rd-order  Runge-Kutta  explicit  scheme  is  used.  The  L, 
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and  L  norms  of  density  error  at  the  solution  points  are  presented  for  tetrahedral  grids  and 
prismatic  grids  in  Table  1  and  2,  respectively.  The  CPR-DG  method  performs  very  well  on  both 
types  of  grid,  achieving  the  nearly  optimal  order  of  accuracy  up  to  6th-order  in  tetrahedral 
meshes  and  4th-order  in  prismatic  meshes. 

5.1.2.  Subsonic  Inviscid  Flow  over  a  Sphere 

In  order  to  verify  the  developed  Euler  solver  on  a  mixed  mesh  with  curved  wall  boundary,  a 
typical  steady  test  case  of  a  subsonic  flow  around  a  sphere  is  considered.  The  freestream  Mach 
number  is  M= 0.3.  Two  computational  grids  are  employed.  One  is  a  purely  prismatic  grid  and  the 
other  is  a  mixed  grid  shown  in  Figures  4  (a)  and  5  (a).  The  mixed  grid  is  composed  of  five  layers 
of  prismatic  cells  around  the  quarter  sphere  and  isotropic  tetrahedral  cells  for  the  remaining 
region.  To  preserve  the  geometry  of  the  sphere  well  with  a  relatively  coarse  mesh,  the  curved 
wall  boundaries  are  represented  by  quadratic  polynomials. 

The  computed  density  contours  obtained  with  the  2nd-  to  4th-order  schemes  are  shown  at 
Figure  4  (b)-(d)  and  Figure  5  (b)-(d).  In  both  grids,  the  trends  of  improvement  in  the  solution  by 
increasing  the  order  of  discretization  are  similar.  The  computed  density  contours  using  the  4th 
order  scheme  appear  to  be  perfectly  symmetric  without  visible  numerical  dissipation  and  also 
quite  smooth  across  the  interface  between  prismatic  and  tetrahedral  cells.  In  this  case,  a  block 
LU-SGS  implicit  scheme  [36,13]  was  used  to  obtain  steady  solutions  efficiently,  and  all  the  cases 
converged  to  machine  zero. 


5.2.  Test  cases  for  the  Navier-Stokes  Equations 
5.2.1.  Accuracy  Study  with  Couette  Flow  Problem 

A  laminar  flow  between  two  parallel  walls  is  considered  here  to  verify  the  discretization  of 
viscous  effects.  The  distance  between  the  walls  is  set  to  H=  10  and  the  computational  domain  is 
chosen  to  be  the  cube  of  [0,  10]x[0,  10]x[0,  10].  The  speed  of  the  moving  upper  wall  (y=1 0)  in 
the  x  direction  is  U=0.3.  The  temperatures  of  the  lower  wall  (v=0)  and  the  upper  one  are  7o=0.8 
and  7]  =0.85  respectively.  The  analytical  solution  for  this  case  is 


u  = — U,  v  =  0,  w  =  0, 

H 

T  =  T0  +  l-(Tl-T0)+^^- 
0  HK  1  °'  2k  H 


1  — 


H 


P=Po .  P  = 


TP 

T' 


(5.2) 


where  yis  specific  heat  ratio  and  k  is  thermal  conductivity.  The  static  pressure  is  set  to  po=Hj 
and  the  viscosity  of  the  fluid  is  assumed  to  be  //=0.01.  The  flow  variables  at  boundary  faces  are 
simply  fixed  to  the  exact  solution. 

Three  successively  refined  prism  grids  are  generated  with  N= 2,  4  and  8  in  the  same  way  as  in 
the  vortex  propagation  case.  Each  cube  is  split  into  two  prisms  by  the  plane  which  is 
perpendicular  to  the  y=0  plane.  The  error  norms  for  the  BR2  formulation  are  presented  in  Table 
3.  The  density  is  used  to  evaluate  the  error.  It  is  shown  that  nearly  the  formal  order  of  accuracy  is 
achieved  for  the  2nd-  to  4th-order  schemes. 
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5.2.2.  Boundary  Layer  Flow  over  a  Flat  Plate 

The  laminar  boundary  layer  flow  over  a  plate  is  then  computed  using  the  CPR  method.  The 
Reynolds  number  based  on  the  plate  length  is  Rex  =10,000  and  the  freestream  Mach  number  is 
M= 0.2.  The  plate  length  L  is  set  to  1.  The  boundary  layer  thickness  at  the  trailing  edge  is 
estimated  by  the  approximate  relation  S=5L/-^Re~x .  The  computational  domain  is  selected  to  be  (- 

2<x<l,  0<y<1005  □  <z<S).  Note  that  the  domain  size  in  the  y-direction  is  chosen  to  be  large 
enough  not  to  affect  the  results  especially  in  the  v- velocity  profiles.  The  freestream  values  are 
specified  at  the  inflow  boundary  at  x=-2  and  the  top  boundary  at  >’=1005  For  the  lower  boundary 
at  y=0,  the  symmetry  conditions  are  used  on  the  upwind  side  to  the  wall  (-2<x<0)  and  the 
adiabatic  wall  conditions  are  imposed  on  the  wall  (0<x<l).  At  the  outflow  boundary  at  x=0,  only 
static  pressure  is  prescribed.  On  the  side  boundaries  at  z=0  and  5,  the  symmetric  conditions  are 
assumed.  First,  we  generated  a  three  dimensional  Cartesian  mesh.  The  grid  cells  are  clustered 
near  the  leading  edge  and  the  cell  sizes  are  increased  geometrically  in  both  x-  and  y-directions.  In 
the  spanwise  z-direction,  we  generate  only  one  cell.  Then  we  divide  each  hexahedral  cell  into 
two  prisms  to  obtain  a  purely  prismatic  grid. 

The  computed  u  and  v  velocity  profiles  are  compared  with  the  Blasius’s  solution  in  Fig.  6. 
The  computational  grid  used  for  the  computations  is  generated  to  have  4  cells  in  the  boundary 
layer  at  x=1.0  and  13  cells  along  the  plate.  The  solution  is  apparently  getting  more  accurate  with 
the  increasing  of  the  order  of  polynomial  approximation,  and  it  is  more  clearly  shown  in  the 
comparison  of  v-profiles.  The  computed  skin  friction  coefficients  on  the  wall  are  also  plotted  at 
Figure  7.  The  agreement  with  the  Blasius’s  solution  also  becomes  better  with  p-order  refinement. 

One  of  the  concerning  issues  when  we  apply  CFD  solver  to  engineering  problems  is  the 
stiffness  arising  from  using  high  aspect  ratio  cells  that  are  clustered  near  the  solid  wall  to  resolve 
the  boundary  layer  especially  in  high  Reynolds  number  flows.  Reynolds  numbers  appearing  in 
aerospace  flow  problems  usually  become  ~106  or  more,  and  so  even  if  we  make  use  of  an 
implicit  time  integration  scheme  for  numerical  simulations,  we  will  likely  encounter  still  small 
time  step  restriction  or  deteriorated  convergence  rate.  A  possible  remedy  for  this  problem  is 
employing  a  line  solver  [27,8].  Here  we  consider  another  approach  to  alleviate  the  stiffness  issue 
by  employing  higher-order  prism  elements  rather  than  having  large  number  of  lower  order 
elements  in  the  boundary  layer.  Since  we  use  a  tensor  basis  polynomial  in  prisms,  we  can  use 
higher  order  polynomial  only  in  the  normal  direction  to  the  wall  while  using  lower  order  one  in 
the  tangential  directions  to  the  wall  so  as  to  prevent  the  unnecessary  increase  of  the 
computational  cost. 

Figure  8  shows  the  computed  Mach  number  by  using  polynomials  of  degree  5  in  the  y- 
direction  and  polynomials  of  degree  2  in  x-  and  z-  directions.  The  grid  has  only  two  cells  in  the 
boundary  layer  at  x=1.0  and  17  cells  along  the  plate.  The  numbers  of  prism  cells  and  DOFs  are 
728  and  26208  respectively.  For  comparison,  we  generated  another  grid  that  has  more  cells  in  the 
boundary  (8  cells  at  x=1.0)  but  the  same  resolution  in  the  x-  and  "-directions  and  employed 
degree  2  polynomials  in  all  directions,  resulting  1736  prisms  and  31248  DOFs.  In  Fig.  9,  the 
computed  v-velocity  profiles  are  shown.  The  computed  profiles  agree  well  with  each  other  and 
also  with  the  Blasius’s  solution.  The  convergence  histories  are  compared  in  Fig.  10.  The 
computations  were  performed  using  the  LU-SGS  scheme  with  the  same  time  step.  Compared  to 
the  computation  using  the  lower  order  scheme  with  the  finer  grid,  employing  the  higher  order 
scheme  with  less  grid  cells  gave  the  reductions  of  about  38%  and  30  %  in  terms  of  Time  steps 
and  CPU  times  to  reach  machine  zero  residual,  although  the  DOFs  are  about  16%  less  than  the 
other’s  one. 
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5.2.3.  Steady  Subsonic  Flow  over  a  Sphere  at  Re=118 

A  steady  viscous  flow  around  a  sphere  is  computed  to  validate  the  developed  NS  solver  on  a 
full  3D  mixed  mesh.  The  Reynolds  number  based  on  the  diameter  was  chosen  to  be  118  so  that 
we  can  compare  the  obtained  results  with  experimental  data  [37]  and  numerical  results  using  the 
SD  scheme  [36,44].  The  Mach  number  is  0.2535  that  is  the  same  value  in  the  reference 
computations.  The  mesh  is  generated  to  have  five  layers  of  prism  cells  and  isotropic  tetrahedral 
cells  for  the  remaining  region.  We  plot  the  cut  of  the  grid  on  a  plane  with  y=0  and  surface  mesh 
on  the  sphere  in  Fig.  11.  The  total  number  of  mixed  cells  is  24,334. 

The  computations  were  performed  using  the  2nd-  to  4th-order  schemes.  The  computed  Mach 
number  contours  and  streamlines  near  wake  using  the  4th-order  CPR  scheme  are  shown  in  Fig. 
12  and  Fig.  13,  respectively.  We  confirmed  that  the  computed  streamlines  and  the  size  of 
separation  region  agree  well  with  both  the  experimental  picture  and  the  numerical  results  in  the 
references.  Here  we  only  show  a  comparison  of  the  computed  skin  friction  profiles  at  the  cross 
section  (y=0)  of  the  sphere  in  Fig.  14.  The  skin  friction  coefficients  computed  by  the  4th-order 
CPR  scheme  and  the  6th-order  SD  scheme  are  right  on  top  of  each  other.  The  3rd-order  CPR 
result  also  agrees  well  with  other  results,  though  one  can  see  only  minor  differences  between 
those  profiles.  The  predicted  separation  angle  using  the  4th-order  CPR  scheme  is  123.6  deg  (the 
wind  side  stagnation  point  has  an  angle  of  0),  which  is  identical  to  the  value  predicted  by  the  6th- 
order  SD  scheme.  In  Fig.  15,  the  computed  drag  coefficient  by  4th-order  CPR  is  compared  to 
available  experimental  data.  The  agreement  is  also  very  good. 

5.2.4.  Unteady  Subsonic  Flow  over  a  Sphere  at  Re=300 

We  consider  an  unsteady  flow  case  over  the  sphere  with  radius  r=l  at  the  Reynolds  number  of 
300  based  on  the  diameter  of  the  sphere.  The  inflow  Mach  number  is  assuemed  to  be  0.3  in  this 
case.  Te  computational  mesh  is  shown  Fig  16.  To  resolve  shedding  vortices,  the  mesh  is 
generated  to  have  finer  cells  in  the  wake  region.  The  total  number  of  mixed  cells  is  54,312.  Local 
grid  size  around  the  sphere  is  -0.2  and  the  size  in  the  wake  region  is  =0.8.  In  this  case,  we 
employed  the  3rd-order  TVD  Runge-Kutta  method  for  the  time  integration  and  computed  by  the 
MPI  parallelized  code  using  8  cores  of  a  cluster  machine  to  reduce  the  wall  clock  time. 

The  computed  Q  isosurface  colored  by  local  Mach  number  using  the  4th-order  CPR  scheme 
is  shown  in  Fig.  17. 

The  obtained  plain  symmetric  wake  vortex  structure  is  comparable  to  the  available  experimental 
and  computational  results  in  [10,20]  at  least  qualitatively.  In  Fig.  18  we  plot  the  history  of  the 
drag  coefficient  Cd  in  terms  of  non-dimensional  time  t.  The  computed  drag  coefficient  and  the 
oscillating  amplitude  of  drag  and  the  Strouhal  number  Str  are  shown  in  the  Table  4.  For 
comparison,  results  from  Gassner  [10]  using  the  4th-order  DG  scheme  on  tetrahedral  grid  and 
from  Tomboulides  [38]  and  Johnson  and  Patel  [20]  obtained  by  incompressible  simulation,  are 
shown  as  well.  The  results  computed  by  the  CPR  method  reasonably  agree  with  those  reference 
values. 
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Table  1.  Test  of  CPR-DG  for  vortex  propagation  problem  (tetrahedral  grids). 


Polynomial 
degree  k 

Grid  size 

L\  error 

L\  order 

Loo  error 

Loo  order 

1 

10x10x10x6 

5.23e-3 

- 

9.56e-2 

- 

20x20x20x6 

1.42e-3 

1.88 

3.57e-2 

1.42 

40x40x40x6 

3.43e-4 

2.05 

9.76e-3 

1.87 

2 

10x10x10x6 

1.68e-3 

- 

6.06e-2 

- 

20x20x20x6 

2.61e-4 

2.69 

1.19e-2 

2.35 

40x40x40x6 

3.77e-5 

2.79 

1.51e-3 

2.98 

3 

10x10x10x6 

4.00e-4 

- 

2.05e-2 

- 

20x20x20x6 

2.44e-5 

4.04 

1.67e-3 

3.62 

40x40x40x6 

1.33e-6 

4.20 

1.00e-4 

4.06 

5 

10x10x10x6 

5.66e-5 

2.34e-3 

20x20x20x6 

9.70e-7 

5.87 

7.78e-5 

4.91 

Table  2.  Test  of  CPR-DG  for  vortex  propagation  problem  (prismatic  grids). 


Polynomial 
degree  k 

Grid 

Li  error 

Li  order 

Loo  error 

Loo  order 

1 

10x10x10x2 

7.37e-3 

- 

1.34e-l 

- 

20x20x20x2 

2.12e-3 

1.80 

4.85e-2 

1.47 

40x40x40x2 

5.19e-4 

2.03 

1.19e-2 

2.03 

2 

10x10x10x2 

2.17e-3 

- 

4.77e-2 

- 

20x20x20x2 

2.67e-4 

3.02 

8.65e-3 

2.46 

40x40x40x2 

2.88e-5 

3.21 

1.04e-3 

3.06 

3 

10x10x10x2 

4.36e-4 

- 

1.54e-2 

- 

20x20x20x2 

2.70e-5 

4.01 

1.43e-3 

3.43 

40x40x40x2 

1.64e-6 

4.04 

9.38e-5 

3.93 
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(a)  Mesh 


(b)  k=  1 
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(c)  k=2 


(d)  k= 3 


Figure  4.  Prismatic  grid  and  computed  density  contours  for  flow  around  a  sphere 


(c)  k= 2 


(d)  k= 3 


Figure  5.  Mixed  grid  (tetrahedrons  and  prisms)  and  computed  density  contours  for  flow  around  a  sphere. 
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Table  3.  Test  of  CPR-DG  (BR2)  for  Couette  flow  problem  (prismatic  grids). 


Polynomial 
degree  k 

Grid 

L\  error 

L\  order 

Loo  error 

Loo  order 

1 

2x2x2x2 

5.5525e-4 

- 

2.4030e-3 

- 

4x4x4x2 

1.1909e-4 

2.221 

3.9968e-4 

2.588 

8x8x8x2 

3.1063e-5 

1.939 

1.1574e-4 

1.788 

2 

2x2x2x2 

8.1732e-6 

- 

2.0928e-5 

- 

4x4x4x2 

1.2867e-6 

2.667 

3.3742e-6 

2.633 

8x8x8x2 

1.6758e-7 

2.941 

5.4916e-7 

2.619 

3 

2x2x2x2 

2.6248e-7 

- 

8.1984e-7 

- 

4x4x4x2 

2.033  le-8 

3.690 

5.7014e-8 

3.846 

8x8x8x2 

1.3907e-9 

3.870 

4.2087e-9 

3.760 

u/U, 


Figure  6.  Comparisons  of  velocity  profiles  in  the  boundary  layer  at  x=0.5. 


u-  and  v-profiles 


on  the  left  and  right. 


Figure  7.  Comparison  the  skin  friction 
coefficient  along  the  plate. 
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Figure  8.  Grid  and  Mach  number  contours  using  the 
CPR  scheme  with  polynomials  of  degree  5  in  the  y- 
direction  (y  direction  stretched  by  factor  10). 


Figure  9.  Comparison  of  v-velocity  profiles  using 
different  degrees  of  polynomial  and  grids. 


Time  step  CPU  [s] 

Figure  10.  Comparisons  of  the  convergence  histories  using  different  degrees  of  polynomial 
and  grids. 
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Figure  11.  Computational  grid  around  a  sphere  for 
the  steady  viscous  flow  over  a  sphere. 


Figure  12.  Computed  pressure  (on  the 
sphere)  and  Mach  number  (on  y=0  plane) 
distributions  using  the  4th-order  CPR 


Figure  13.  Computed  streamlines  using 
the  4th-order  CPR  scheme  near  the  wake 
region  behind  the  sphere. 
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Figure  14.  Comparison  of  computed  skin 
friction  coefficients  using  the  3rd-  and  4th- 
order  CPR  schemes  and  6th-order  SD 


Re 

Figure  15.  Comparison  between  the 
computed  drag  coefficient  using  the  4th- 
order  CPR  scheme  and  experimental  data 


Figure  16.  Computational  grid  around  a  sphere  for  the  unsteady  viscous  flow  over  a  sphere. 
(Left:  Entire  grid,  Right:  Grid  around  the  sphere.) 
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Figure  17.  Computed  Q  isosurface  in  the 


wake  region  of  the  viscous  laminar  flow 
over  a  sphere  at  Re=300. 
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Figure  18.  Time  history  of  the  drag 
coefficient. 


Table  4.  Comparisons  of  the  averaged  drag  coefficient,  the  amplitude  of  drag  and  the  Strouhal 

number. 


Method 

Cd 

ACd 

Str 

Present  (LCP) 

0.670 

0.0032 

0.131 

Gassner  [10] 

0.673 

0.0031 

0.135 

Tomboulides  [38] 

0.671 

0.0028 

0.136 

Johnson  &  Patel 
[20] 

0.656 

0.0035 

0.137 

6.  Conclusion 

The  CPR  method  is  successfully  extended  to  3D  hybrid  unstructured  meshes  using  tetrahedral 
and  prismatic  elements.  The  CPR  formulation  for  tetrahedral  elements  is  directly  derived  in  the 
same  manner  for  2D  triangular  elements  and  the  one  for  prism  is  obtained  by  just  a  combination 
of  the  ID  and  2D  schemes.  The  resulting  scheme  needs  no  explicit  integrations  and  no  data 
reconstructions.  This  numerical  efficiency  is  more  significant  in  3D  simulations  in  comparison  to 
2D  simulations  because  numerical  complexities  involved  in  high-order  quadratures  and 
reconstructions  rapidly  increase  in  3D. 

The  developed  CPR  scheme  is  verified  with  grid  convergence  studies  for  an  inviscid  flow  and 
a  viscous  flow,  indicating  that  the  developed  scheme  is  capable  of  achieving  nearly  the  optimal 
order  of  accuracy.  Then,  several  validation  cases  are  computed  for  solving  the  3D  Euler 
equations  and  the  3D  NS  equations.  The  CPR  method  performs  very  well  to  obtain  high-order 
accurate  solution  for  all  cases.  Future  studies  include  extension  to  adopt  hexahedral  and 
pyramidal  cells  for  more  flexible  geometry  discretizations  and  hp-adaptation  techniques  for 
realizing  practical  high  accurate  CFD  simulations. 
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